home *** CD-ROM | disk | FTP | other *** search
/ Software Vault: The Diamond Collection / The Diamond Collection (Software Vault)(Digital Impact).ISO / cdr44 / newmat08.zip / NEWMAT2.CPP < prev    next >
C/C++ Source or Header  |  1995-01-15  |  14KB  |  439 lines
  1. //$$ newmat2.cpp      Matrix row and column operations
  2.  
  3. // Copyright (C) 1991,2,3,4: R B Davies
  4.  
  5. #define WANT_MATH
  6.  
  7. #include "include.h"
  8.  
  9. #include "newmat.h"
  10. #include "newmatrc.h"
  11.  
  12. //#define REPORT { static ExeCounter ExeCount(__LINE__,2); ++ExeCount; }
  13.  
  14. #define REPORT {}
  15.  
  16. //#define MONITOR(what,storage,store) \
  17. //   { cout << what << " " << storage << " at " << (long)store << "\n"; }
  18.  
  19. #define MONITOR(what,store,storage) {}
  20.  
  21. /************************** Matrix Row/Col functions ************************/
  22.  
  23. void MatrixRowCol::Add(const MatrixRowCol& mrc)
  24. {
  25.    REPORT
  26.    int f = mrc.skip; int l = f + mrc.storage; int lx = skip + storage;
  27.    if (f < skip) f = skip; if (l > lx) l = lx; l -= f;
  28.    if (l<=0) return;
  29.    Real* elx=store+f; Real* el=mrc.store+f;
  30.    while (l--) *elx++ += *el++;
  31. }
  32.  
  33. void MatrixRowCol::AddScaled(const MatrixRowCol& mrc, Real x)
  34. {
  35.    REPORT
  36.    int f = mrc.skip; int l = f + mrc.storage; int lx = skip + storage;
  37.    if (f < skip) f = skip; if (l > lx) l = lx; l -= f;
  38.    if (l<=0) return;
  39.    Real* elx=store+f; Real* el=mrc.store+f;
  40.    while (l--) *elx++ += *el++ * x;
  41. }
  42.  
  43. void MatrixRowCol::Sub(const MatrixRowCol& mrc)
  44. {
  45.    REPORT
  46.    int f = mrc.skip; int l = f + mrc.storage; int lx = skip + storage;
  47.    if (f < skip) f = skip; if (l > lx) l = lx; l -= f;
  48.    if (l<=0) return;
  49.    Real* elx=store+f; Real* el=mrc.store+f;
  50.    while (l--) *elx++ -= *el++;
  51. }
  52.  
  53. void MatrixRowCol::Inject(const MatrixRowCol& mrc)
  54. // copy stored elements only
  55. {
  56.    REPORT
  57.    int f = mrc.skip; int l = f + mrc.storage; int lx = skip + storage;
  58.    if (f < skip) f = skip; if (l > lx) l = lx; l -= f;
  59.    if (l<=0) return;
  60.    Real* elx = store+f; Real* ely = mrc.store+f;
  61.    while (l--) *elx++ = *ely++;
  62. }
  63.  
  64. Real DotProd(const MatrixRowCol& mrc1, const MatrixRowCol& mrc2)
  65. {
  66.    REPORT                                         // not accessed
  67.    int f = mrc1.skip; int f2 = mrc2.skip;
  68.    int l = f + mrc1.storage; int l2 = f2 + mrc2.storage;
  69.    if (f < f2) f = f2; if (l > l2) l = l2; l -= f;
  70.    if (l<=0) return 0.0;
  71.    
  72.    Real* el1=mrc1.store+f; Real* el2=mrc2.store+f;
  73.    Real sum = 0.0;
  74.    while (l--) sum += *el1++ * *el2++;
  75.    return sum;
  76. }
  77.  
  78. void MatrixRowCol::Add(const MatrixRowCol& mrc1, const MatrixRowCol& mrc2)
  79. {
  80.    int f = skip; int l = skip + storage;
  81.    int f1 = mrc1.skip; int l1 = f1 + mrc1.storage;
  82.    if (f1<f) f1=f; if (l1>l) l1=l;
  83.    int f2 = mrc2.skip; int l2 = f2 + mrc2.storage;
  84.    if (f2<f) f2=f; if (l2>l) l2=l;
  85.    Real* el = store + f;
  86.    Real* el1 = mrc1.store+f1; Real* el2 = mrc2.store+f2;
  87.    if (f1<f2)
  88.    {
  89.       int i = f1-f; while (i--) *el++ = 0.0;
  90.       if (l1<=f2)                              // disjoint
  91.       {
  92.      REPORT                                // not accessed
  93.          i = l1-f1;     while (i--) *el++ = *el1++;
  94.          i = f2-l1;     while (i--) *el++ = 0.0;
  95.          i = l2-f2;     while (i--) *el++ = *el2++;
  96.          i = l-l2;      while (i--) *el++ = 0.0;
  97.       }
  98.       else
  99.       {
  100.          i = f2-f1;    while (i--) *el++ = *el1++;
  101.          if (l1<=l2)
  102.          {
  103.             REPORT
  104.             i = l1-f2; while (i--) *el++ = *el1++ + *el2++;
  105.             i = l2-l1; while (i--) *el++ = *el2++;
  106.             i = l-l2;  while (i--) *el++ = 0.0;
  107.          }
  108.          else
  109.          {
  110.             REPORT
  111.             i = l2-f2; while (i--) *el++ = *el1++ + *el2++;
  112.             i = l1-l2; while (i--) *el++ = *el1++;
  113.             i = l-l1;  while (i--) *el++ = 0.0;
  114.          }
  115.       }
  116.    }
  117.    else
  118.    {
  119.       int i = f2-f; while (i--) *el++ = 0.0;
  120.       if (l2<=f1)                              // disjoint
  121.       {
  122.      REPORT                                // not accessed
  123.          i = l2-f2;     while (i--) *el++ = *el2++;
  124.          i = f1-l2;     while (i--) *el++ = 0.0;
  125.      i = l1-f1;     while (i--) *el++ = *el1++;
  126.      i = l-l1;      while (i--) *el++ = 0.0;
  127.       }
  128.       else
  129.       {
  130.          i = f1-f2;    while (i--) *el++ = *el2++;
  131.          if (l2<=l1)
  132.          {
  133.         REPORT
  134.             i = l2-f1; while (i--) *el++ = *el1++ + *el2++;
  135.             i = l1-l2; while (i--) *el++ = *el1++;
  136.             i = l-l1;  while (i--) *el++ = 0.0;
  137.          }
  138.          else
  139.          {
  140.         REPORT
  141.             i = l1-f1; while (i--) *el++ = *el1++ + *el2++;
  142.             i = l2-l1; while (i--) *el++ = *el2++;
  143.             i = l-l2;  while (i--) *el++ = 0.0;
  144.          }
  145.       }
  146.    }
  147. }
  148.  
  149. void MatrixRowCol::Sub(const MatrixRowCol& mrc1, const MatrixRowCol& mrc2)
  150. {
  151.    int f = skip; int l = skip + storage;
  152.    int f1 = mrc1.skip; int l1 = f1 + mrc1.storage;
  153.    if (f1<f) f1=f; if (l1>l) l1=l;
  154.    int f2 = mrc2.skip; int l2 = f2 + mrc2.storage;
  155.    if (f2<f) f2=f; if (l2>l) l2=l;
  156.    Real* el = store + f;
  157.    Real* el1 = mrc1.store+f1; Real* el2 = mrc2.store+f2;
  158.    if (f1<f2)
  159.    {
  160.       int i = f1-f; while (i--) *el++ = 0.0;
  161.       if (l1<=f2)                              // disjoint
  162.       {
  163.      REPORT                                // not accessed
  164.          i = l1-f1;     while (i--) *el++ = *el1++;
  165.          i = f2-l1;     while (i--) *el++ = 0.0;
  166.          i = l2-f2;     while (i--) *el++ = - *el2++;
  167.          i = l-l2;      while (i--) *el++ = 0.0;
  168.       }
  169.       else
  170.       {
  171.          i = f2-f1;    while (i--) *el++ = *el1++;
  172.          if (l1<=l2)
  173.          {
  174.         REPORT
  175.             i = l1-f2; while (i--) *el++ = *el1++ - *el2++;
  176.             i = l2-l1; while (i--) *el++ = - *el2++;
  177.             i = l-l2;  while (i--) *el++ = 0.0;
  178.          }
  179.          else
  180.          {
  181.         REPORT
  182.             i = l2-f2; while (i--) *el++ = *el1++ - *el2++;
  183.             i = l1-l2; while (i--) *el++ = *el1++;
  184.             i = l-l1;  while (i--) *el++ = 0.0;
  185.          }
  186.       }
  187.    }
  188.    else
  189.    {
  190.       int i = f2-f; while (i--) *el++ = 0.0;
  191.       if (l2<=f1)                              // disjoint
  192.       {
  193.      REPORT                                // not accessed
  194.          i = l2-f2;     while (i--) *el++ = - *el2++;
  195.          i = f1-l2;     while (i--) *el++ = 0.0;
  196.          i = l1-f1;     while (i--) *el++ = *el1++;
  197.          i = l-l1;      while (i--) *el++ = 0.0;
  198.       }
  199.       else
  200.       {
  201.          i = f1-f2;    while (i--) *el++ = - *el2++;
  202.          if (l2<=l1)
  203.          {
  204.         REPORT
  205.             i = l2-f1; while (i--) *el++ = *el1++ - *el2++;
  206.             i = l1-l2; while (i--) *el++ = *el1++;
  207.             i = l-l1;  while (i--) *el++ = 0.0;
  208.          }
  209.          else
  210.          {
  211.             REPORT
  212.             i = l1-f1; while (i--) *el++ = *el1++ - *el2++;
  213.             i = l2-l1; while (i--) *el++ = - *el2++;
  214.             i = l-l2;  while (i--) *el++ = 0.0;
  215.          }
  216.       }
  217.    }
  218. }
  219.  
  220.  
  221. void MatrixRowCol::Add(const MatrixRowCol& mrc1, Real x)
  222. {
  223.    REPORT
  224.    int f = mrc1.skip; int l = f + mrc1.storage; int lx = skip + storage;
  225.    if (f < skip) { f = skip; if (l < f) l = f; }
  226.    if (l > lx) { l = lx; if (f > lx) f = lx; }
  227.  
  228.    Real* elx = store+skip; Real* ely = mrc1.store+f;
  229.  
  230.    int l1 = f-skip;  while (l1--) *elx++ = x;
  231.        l1 = l-f;     while (l1--) *elx++ = *ely++ + x;
  232.        lx -= l;      while (lx--) *elx++ = x;
  233. }
  234.  
  235. void MatrixRowCol::RevSub(const MatrixRowCol& mrc1)
  236. {
  237.    REPORT
  238.    int f = mrc1.skip; int l = f + mrc1.storage; int lx = skip + storage;
  239.    if (f < skip) { f = skip; if (l < f) l = f; }
  240.    if (l > lx) { l = lx; if (f > lx) f = lx; }
  241.  
  242.    Real* elx = store+skip; Real* ely = mrc1.store+f;
  243.  
  244.    int l1 = f-skip;  while (l1--) { *elx = - *elx; elx++; }
  245.        l1 = l-f;     while (l1--) { *elx = *ely++ - *elx; elx++; }
  246.        lx -= l;      while (lx--) { *elx = - *elx; elx++; }
  247. }
  248.  
  249. void MatrixRowCol::ConCat(const MatrixRowCol& mrc1, const MatrixRowCol& mrc2)
  250. {
  251.    REPORT
  252.    int f1 = mrc1.skip; int l1 = f1 + mrc1.storage; int lx = skip + storage;
  253.    if (f1 < skip) { f1 = skip; if (l1 < f1) l1 = f1; }
  254.    if (l1 > lx) { l1 = lx; if (f1 > lx) f1 = lx; }
  255.  
  256.    Real* elx = store+skip;
  257.  
  258.    int i = f1-skip;  while (i--) *elx++ =0.0;
  259.    i = l1-f1;
  260.    if (i)                       // in case f1 would take ely out of range
  261.       { Real* ely = mrc1.store+f1;  while (i--) *elx++ = *ely++; }
  262.  
  263.    int f2 = mrc2.skip; int l2 = f2 + mrc2.storage; i = mrc1.length;
  264.    int skipx = l1 - i; lx -= i; // addresses rel to second seg, maybe -ve
  265.    if (f2 < skipx) { f2 = skipx; if (l2 < f2) l2 = f2; }
  266.    if (l2 > lx) { l2 = lx; if (f2 > lx) f2 = lx; }
  267.  
  268.    i = f2-skipx; while (i--) *elx++ = 0.0;
  269.    i = l2-f2;
  270.    if (i)                       // in case f2 would take ely out of range
  271.       { Real* ely = mrc2.store+f2; while (i--) *elx++ = *ely++; }
  272.    lx -= l2;     while (lx--) *elx++ = 0.0;
  273. }
  274.  
  275. void MatrixRowCol::Multiply(const MatrixRowCol& mrc1)
  276. // element by element multiply into
  277. {
  278.    REPORT
  279.    int f = mrc1.skip; int l = f + mrc1.storage; int lx = skip + storage;
  280.    if (f < skip) { f = skip; if (l < f) l = f; }
  281.    if (l > lx) { l = lx; if (f > lx) f = lx; }
  282.  
  283.    Real* elx = store+skip; Real* ely = mrc1.store+f;
  284.  
  285.    int l1 = f-skip;  while (l1--) *elx++ = 0;
  286.        l1 = l-f;     while (l1--) *elx++ *= *ely++;
  287.        lx -= l;      while (lx--) *elx++ = 0;
  288. }
  289.  
  290. void MatrixRowCol::Multiply(const MatrixRowCol& mrc1, const MatrixRowCol& mrc2)
  291. // element by element multiply
  292. {
  293.    int f = skip; int l = skip + storage;
  294.    int f1 = mrc1.skip; int l1 = f1 + mrc1.storage;
  295.    if (f1<f) f1=f; if (l1>l) l1=l;
  296.    int f2 = mrc2.skip; int l2 = f2 + mrc2.storage;
  297.    if (f2<f) f2=f; if (l2>l) l2=l;
  298.    Real* el = store + f; int i;
  299.    if (f1<f2) f1 = f2; if (l1>l2) l1 = l2;
  300.    if (l1<=f1) { REPORT i = l-f; while (i--) *el++ = 0.0; }  // disjoint
  301.    else
  302.    {
  303.       REPORT
  304.       Real* el1 = mrc1.store+f1; Real* el2 = mrc2.store+f1;
  305.       i = f1-f ;    while (i--) *el++ = 0.0;
  306.       i = l1-f1;    while (i--) *el++ = *el1++ * *el2++;
  307.       i = l-l1;     while (i--) *el++ = 0.0;
  308.    }
  309. }
  310.  
  311. void MatrixRowCol::Copy(const MatrixRowCol& mrc1)
  312. {
  313.    REPORT
  314.    int f = mrc1.skip; int l = f + mrc1.storage; int lx = skip + storage;
  315.    if (f < skip) { f = skip; if (l < f) l = f; }
  316.    if (l > lx) { l = lx; if (f > lx) f = lx; }
  317.  
  318.    Real* elx = store+skip; Real* ely = mrc1.store+f;
  319.  
  320.    int l1 = f-skip;  while (l1--) *elx++ = 0.0;
  321.        l1 = l-f;     while (l1--) *elx++ = *ely++;
  322.        lx -= l;      while (lx--) *elx++ = 0.0;
  323. }
  324.  
  325. void MatrixRowCol::CopyCheck(const MatrixRowCol& mrc1)
  326. // Throw an exception this would lead to a loss of data
  327. {
  328.    REPORT
  329.    int f = mrc1.skip; int l = f + mrc1.storage; int lx = skip + storage;
  330.    if (f < skip || l > lx) Throw(ProgramException("Illegal Conversion"));
  331.  
  332.    Real* elx = store+skip; Real* ely = mrc1.store+f;
  333.  
  334.    int l1 = f-skip;  while (l1--) *elx++ = 0.0;
  335.        l1 = l-f;     while (l1--) *elx++ = *ely++;
  336.        lx -= l;      while (lx--) *elx++ = 0.0;
  337. }
  338.  
  339. void MatrixRowCol::Negate(const MatrixRowCol& mrc1)
  340. {
  341.    REPORT
  342.    int f = mrc1.skip; int l = f + mrc1.storage; int lx = skip + storage;
  343.    if (f < skip) { f = skip; if (l < f) l = f; }
  344.    if (l > lx) { l = lx; if (f > lx) f = lx; }
  345.  
  346.    Real* elx = store+skip; Real* ely = mrc1.store+f;
  347.  
  348.    int l1 = f-skip;  while (l1--) *elx++ = 0.0;
  349.        l1 = l-f;     while (l1--) *elx++ = - *ely++;
  350.        lx -= l;      while (lx--) *elx++ = 0.0;
  351. }
  352.  
  353. void MatrixRowCol::Multiply(const MatrixRowCol& mrc1, Real s)
  354. {
  355.    REPORT
  356.    int f = mrc1.skip; int l = f + mrc1.storage; int lx = skip + storage;
  357.    if (f < skip) { f = skip; if (l < f) l = f; }
  358.    if (l > lx) { l = lx; if (f > lx) f = lx; }
  359.  
  360.    Real* elx = store+skip; Real* ely = mrc1.store+f;
  361.  
  362.    int l1 = f-skip;  while (l1--) *elx++ = 0.0;
  363.        l1 = l-f;     while (l1--) *elx++ = *ely++ * s;
  364.        lx -= l;      while (lx--) *elx++ = 0.0;
  365. }
  366.  
  367. void DiagonalMatrix::Solver(MatrixRowCol& mrc, const MatrixRowCol& mrc1)
  368. {
  369.    REPORT
  370.    int f = mrc1.skip; int f0 = mrc.skip;
  371.    int l = f + mrc1.storage; int lx = f0 + mrc.storage;
  372.    if (f < f0) { f = f0; if (l < f) l = f; }
  373.    if (l > lx) { l = lx; if (f > lx) f = lx; }
  374.  
  375.    Real* elx = mrc.store+f0; Real* eld = store+f;
  376.  
  377.    int l1 = f-f0;    while (l1--) *elx++ = 0.0;
  378.        l1 = l-f;     while (l1--) *elx++ /= *eld++;
  379.        lx -= l;      while (lx--) *elx++ = 0.0;
  380.    // Solver makes sure input and output point to same memory
  381. }
  382.  
  383. void MatrixRowCol::Copy(const Real*& r)
  384. {
  385.    REPORT
  386.    Real* elx = store+skip; const Real* ely = r+skip; r += length;
  387.    int l = storage; while (l--) *elx++ = *ely++;
  388. }
  389.  
  390. void MatrixRowCol::Copy(Real r)
  391. {
  392.    REPORT
  393.    Real* elx = store+skip; int l = storage; while (l--) *elx++ = r;
  394. }
  395.  
  396. Real MatrixRowCol::SumAbsoluteValue()
  397. {
  398.    REPORT
  399.    Real sum = 0.0; Real* elx = store+skip; int l = storage;
  400.    while (l--) sum += fabs(*elx++);
  401.    return sum;
  402. }
  403.  
  404. Real MatrixRowCol::Sum()
  405. {
  406.    REPORT
  407.    Real sum = 0.0; Real* elx = store+skip; int l = storage;
  408.    while (l--) sum += *elx++;
  409.    return sum;
  410. }
  411.  
  412. void MatrixRowCol::SubRowCol(MatrixRowCol& mrc, int skip1, int l1) const
  413. {
  414.    mrc.length = l1;  mrc.store = store + skip1;  int d = skip - skip1;
  415.    mrc.skip = (d < 0) ? 0 : d;  d = skip + storage - skip1;
  416.    d = ((l1 < d) ? l1 : d) - mrc.skip;  mrc.storage = (d < 0) ? 0 : d;
  417.    mrc.cw = 0;
  418. }
  419.  
  420. MatrixRowCol::~MatrixRowCol()
  421. {
  422.    int t1 = +(cw*IsACopy); int t2 = !(cw*StoreHere);
  423.    if (t1 && t2)
  424.    {
  425.       Real* f = store+skip;
  426.       MONITOR_REAL_DELETE("Free    (RowCol)",storage,f) 
  427. #ifdef Version21
  428.       delete [] f;
  429. #else
  430.       delete [storage] f;
  431. #endif
  432.    }
  433. }
  434.  
  435. MatrixRow::~MatrixRow() { if (+(cw*StoreOnExit)) gm->RestoreRow(*this); }
  436.  
  437. MatrixCol::~MatrixCol() { if (+(cw*StoreOnExit)) gm->RestoreCol(*this); }
  438.  
  439.